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ABSTRACT 

The  eigenvector  solution  of  the  time-invariant  matrix  Riccati  equation 
is  discussed.   The  coefficient  matrix  of  the  canonical  equation  is  allowed 
to  have  multiple  eigenvalues,  namely,  the  matrix  could  be  either  derog- 
atory or  defective.   The  solution  of  matrix  Riccati  equation  is  then  cal- 
culated from  a  part  of  similarity  transformation  which  should  reduce 
the  coefficient  matrix  to  the  Jordan  canonical  form. 
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1 .   INTRODUCTION 

It  is  well  known  [1,2]  that  the  linear  regulator  problem  with  quadratic 
cost  functional  is  reduced  to  the  problem  of  the  matrix  Riccati  equation. 

Although  the  solution  of  the  matrix  Riccati  equation  exists  and  is 
unique  [2,3],  it  is  not  easy  to  obtain  a  numerical  solution  for  the 
high-order  system. 

One  of  the  numerical  methods  is  the  eigenvector  solution  which  has  been 
studied  by  several  authors  [U , 5 »6 ] .   In  their  studies,  however,  the  coeffi- 
cient matrix  of  the  canonical  equation  is  assumed  to  have  distinct  eigen- 
values, which  is  a  restriction  to  the  method  of  eigenvector  solutions. 
Recently,  Martensson  [7]  discussed  the  case  in  which  the  coefficient  matrix 
is  nondiagonable.   Similarly,  in  this  paper,  the  assumption  of  distinct 
eigenvalues  will  also  be  removed  so  that  the  coefficient  matrix  is  allowed 
to  have  the  multiple  eigenvalues,  namely,  the  matrix  is  either  derogatory 
or  defective.   To  do  this,  we  explicitly  construct  the  similarity  transfor- 
mation which  should  reduce  the  coefficient  matrix  to  the  Jordan  canonical 
form,  either  diagonal  or  nondiagonal,  and  then  the  solution  of  the  matrix 
Riccati  equation  is  calculated  from  a  part  of  the  similarity  transformation 
matrix. 

Since  matrix  Riccati  equations  of  fairly  large  sizes  appear  in  many 
engineering  applications  and  since  the  solution  algorithm  can  be  easily 
constructed  such  that  it  becomes  suitable  for  a  parallel  computer  [8],  we 
provide  in  Appendix  II  the  basic  codes  written  in  ILLIAC  IV  language, 
GLYPNIR. 


2.   FORMULATION  OF  THE  PROBLEM 

In  this  section,  the  source  of  the  matrix  Riccati  equation  and  its 
relation  to  optimal  control  problems  are  summarized  in  short  for  our 
further  discussions  [1,2,3]. 

Let  us  consider  the  linear  time- invariant  dynamical  system 

x(t)  =  Ax(t)  +  Bu(t) 

y(t)  =  Cx(t)  (2.1) 

x(tQ)=x0 
where  A,  B,  and  C  are  n  x  n,  n  x  r,  and  m  x  n  constant  matrices,  respec- 
tively, and  x(t)  is  an  n-dimensional  state  vector,  u(t)  is  an  r-dimensional 
control  vector,  and  y(t)  is  an  m-dimensional  output  vector,  respectively. 
Further,  we  assume  that  the  system  (2.1)  is  completely  observable.   The 
optimal  output  regulator  problem  is  then  to  determine  the  control  u(t)  which 
minimizes  the  quadratic  cost  functional, 

V  =  |yT  (t  )   IV(t  )  +  \       Jf  (yT(t)  Qy(t)  +  uT(t)Ru(t))dt 

t0  (2.2) 

where  P  and  Q  are  m  x  m  symmetric  positive   semidefinite  matrices   and  R  is 
an  r  x  r  symmetric  positive   definite  matrix.      Substituting  y(t)   =   Cx(t) 
into    (2.2),    V     is   rewritten   as 

V  =  \  xT(tf )    CTPCx(tf)    +  |     /      (xT(t)CTQC  x(t)    +  uT(t)Ru(t))dt 

t0  (2.3) 

T        T 
Since  the  system  (2.1)  is  completely  observable,  C  PC  and  C  QC  are  symmetric 

and  positive  semidefinite  [l]. 

The  optimal  feedback  control,  therefore,  is  given  by 

u*(t)  =  -R_1BTK(t)x(t)  -  (2.U) 

where  K(t)  is  an  n  x  n  symmetric  positive  definite  matrix  which  is  the 


solution  of  the  matrix  Riccati  equation 

K(t)  +  K(t)  A  +  ATK(t)  -  K(t)BR_1BTK(t)  +  CTQC  =  0         (2.5) 
with  the  boundary  condition 

K(t  )  =  CTPC  '  (2.6) 

Furthermore,  the  optimal  trajectory  is  the  solution  of  the  system  of 
differential  equations 

x(t)  =  (A  -  BR"1BTK(t))  x(t)  (2.7) 

x(tQ)  =  xQ 

The  minimum  cost  is  given  by 

V  •  =  |xT(t)  K(t)  x(t)  (2.8) 


In  addition  to  the  assumption  of  complete  observability  we  further 
assume  that  the  system  is  completely  controllable.   Setting  P  =  0  and 
t  =  °°,  the  output  regulator  problem  turns  out  to  be  the  following: 
minimize 

V  =  \     /~(yT(t)Qy(t)  +  uT(t)Ru(t))  dt  (2.9) 

subject  to 

x(t)  =  Ax(t)  +  Bu(t) 

y(t)  =  Cx(t)  (2.10) 

x(0)  =  xQ 
Then  the  optimal  control  is  given  by 

u(t)  =  -R_1BTK  x(t)  (2.11) 

where  K  is  the  n  x  n  symmetric  positive  definite  matrix  which  is  the  solu- 
tion of  algebraic  Riccati  equation 

KA  +  ATK  -  KBR_1BTK  +  CTQC  =  0  (2.12) 

The  optimal  trajectory  is  the  solution  of  the  system 


x(t)  =  Gx(t) 

-1  T 
G  =  A  -  BR  B  K 


(2.13) 


x 


(0)  = 


Note  that  for  an  optimal  control  system  the  matrix  G  is  diagonable  and  has' 
eigenvalues  with  negative  real  parts.   Therefore,  the  symmetric  positive 
definite  solution  K  of  (2.12)  must  be  such'  that  G  has  eigenvalues  with 
negative  real  parts. 

Kalman  [2]  has  shown  that  if  the  system  (2.1)  is  completely  observable 
and  completely  controllable,  then 

lim  K(t  )  =  K  (2.1U) 

t  -*» 

f 

We   shall  observe  this   equation  in  the   following. 

Let   K(t)   =  Z(t)X-1(t)  (2.15) 

where  X(t)    and  Z(t)    are  n  x  n  matrices    and  X(t)    is   non-singular.      Substituting 

K(t)   =    (Z(t)   -  K(t)   X(t))   X_1(t)  (2.16) 

into    (2.5)    and  setting 

X(t)   =  AX(t)    -   BR-1BT   Z(t)  (2.17) 

we  obtain 

Z(t)    =   -CTQCX(t)   -  ATZ(t)  (2.18) 

Therefore    (2.5)  becomes    a  system  of  2n-dimensional   linear  homogeneous 
differential   equations 

x(t) 


"i(t)  1=  \\ 

z(t)      [-c1 


-1  T 
-BR  B 


-C  QC    -A 


Z(t) 


(2.19) 


Since  we  are  interested  in  increasing  the  time  t,  let 


t  =  tf  -  t 


(2.20) 


Then 


where 


X(r) 

=  ¥ 

"*(*) 

_Z(t)_ 

Z(t) 

-A 

-it" 

BR      B 

¥  = 

T 
C   QC 

AT 

„ 

—* 

(2.21) 


with 

X(0)  =  I 

Z(0)  =  CTQC  (2.22) 

Therefore,  the  solution  will  be  given  by  the  following  matrix  exponential 

X(f)" 


Z(t) 


=  e 


¥t 


T 
CTQC 


(2.23) 


Consequently,  the  solution  of  the  matrix  Riccati  equation  will  be  obtained 
by  (2.15)  as  t->  °°.   In  the  next  section,  we  shall  show  that  two  matrices 
X(t)  and  Z(t)  are  obtained  from  the  similarity  transformation  which  reduces 
¥  into  the  Jordan  canonical  form.   Throughout . the  remaining  discussions, 
the  increasing  time  variable  t  will  be  used  instead  of  t. 


3.   SOLUTION  OF  THE  MATRIX  RICCATI  EQUATION 


We  begin  with  two  lemmas  for  our  main  result. 


Lemma  1.   The  eigenvalues  of  the  matrix  ¥  of  (2.21)  are  symmetric 
with  respect  to  the  imaginary  axis. 


Proof.   Let  I  be  a  2n  x  2n  matrix  given  by 


xi  = 


-I 
0 


(3.1) 


where   I   is   an  n  x  n   identity  matrix.      Then   it    is  obvious  that 

-1  T 

I  =   I 

1  1 

T  T 

I  WI        =   -¥ 


(3.2) 
(3.3) 


Since  the  eigenvalues  of  a  matrix  are  unchanged  by  either  similarity  trans- 
formations or  transposing,  the  eigenvalues  of  W  occur  in  pairs  with  opposite 


signs. 


Lemma  2.   Let 


T  = 


All   A12 


A      A 
LA21   A22J 


(3.10 


where  A   ,  A   ,  A   ,  and  A   are  n  x  n,  m  x  m,  nxm,  and  m  x  n  matrices, 
respectively.   Assume  that  A.^   is  non-singular.   Then 


-1 


T  T 

11  12 


T  T 

L     21  22J 

-1 


(3.5) 


exists    iff  A        -  A        A        "  A        is   non-singular,    and  is   given  by 


-1 


-1 


Tll  :  All  "  +  All  "  ^2  T22  A21  All 
T12  =  "All  '  A12  T22 


-1 


(3.6) 


T    =  -T    A    A 
21     22   21   11 


-1 


T22  =  (A22  "  A21  All    A12) 


7 


where  submatrices  T   ,  T   ,  T   ,  and  T   are  the  same  sizes  as  A   , 
A22'  A12'  and  A21'  resPectively- 

Proof.   Suppose  A   -  A   Ann"  ^o  is  non-sin8'ular'   From  TT  "  =  I. 

All  Tll  +  A12  T21  =  Xn 


(3.7) 


lll  T12  +  ^2  T22  "  ° 


A21  Tll  +  A22  T21  "  ° 

A21  T12  +  A22  T22  =  Im 
Premultiplying  the  second  equation  by  A   A     and  subtracting  from  the 


fourth, 


T22  "  (A22  "  A21  All  "  A12)  ' 


And  hence 


T12  =   ^1   \2   T22 


Similarly,  from  the  first  and  third, 
T21  =  "T22  A21  ^1  ' 


-1 


-1 


Tll  ==  All  "  +  All  "  h2   T22  A21  All 


-1 


Hence  T        is   a  right    inverse  of  T.      In   a  similar  way  we   can  show  that   T 
is   also   a  left    inverse  of  T.      Conversely,    suppose  T   is   non-singular.      Then 


0  ?     det 


A, 


^1     "12 


LA21     A22J 


=  det 


A22~A21  hi        A12. 


-1 


-det    (A^)    det(A22-A21A11-     A^) 


(3.8) 


Therefore,  A   -  A   A     A   is  nonsingular. 


Assuming  that  W  has  distinct  eigenvalues,  we  now  turn  to  constructing 
the  similarity  transformation  with  which  W  of  (2.2l)  is  reduced  to  the 


8 


diagonal  form.   By  lemma  1,  W  has  n  eigenvalues  with  positive  real  parts 
and  n  eigenvalues  with  negative  real  parts.   If  we  construct  the  similarity- 
transformation  S,  whose  columns  are  the  right  eigenvectors  of  ¥,  such  that 
the  first  n  columns  are  the  right  eigenvectors  of  W  corresponding  to  the 
n  eigenvalues  of  W  with  positive  real  parts,  then 


s_1ws  = 


1 


X. 


*A 


2n 


(3.9) 


,-1 


where,    of  course,   the  rows   of  S  "   are  the  left   eigenvectors   of  W. 
O'donnell    [5]   and  Potter    [6]  have   shown  that   if  we  partition  S   into   four 
n  x  n  submatrices, 


S  = 


Sll        S12 


.321        S22J 


(3.10) 


then  the  solution  of  the  matrix  Riccati  equation  is  given  by 

-1 


K  -   S21   S±1 
provided  that   S        is  non-singular. 


(3.11) 


We   shall  now  remove  the   restriction  that  W  has   distinct   eigenvalues   and 
present   an  alternative  proof  to  that   of  Martensson    [7]  to   show  that    (3.11 ) 
holds    for  the   case   in  which  ¥  has   nondiagonal  Jordan   form. 

Let   us    consider  the   general  Jordan   form  J  with  the   similarity  transfor- 
mation S,   which  will  be   constructed  later, 
S-1¥S  =  J 


where,      J  = 


"i , 


(3.12) 


'J 


m 


and  in  which 


m. 


A.    1 

1 


A.    1 

l 


".  A. 


,  (m.  x  m.  ) 


(3.13) 


If  we  denote  S  by  its  column  vectors  s  ,  . . . ,  s   ,  then 


¥  Ls1,  s2, 


>    S2n]  =  [S1'  S2'  ••••  S2n]  J 


(3.1M 


For  any  Jordan  block  J   ,  if  i.  >  1  then,  dropping  subscripts  of  s  for 

m.      l 

l 

convenience, 

(1)   .     _(D 


Ws 


=  sv    'A. 

l 


w«<2>    ■   *(l)   +  >      .(2) 
ws  =   s  +   A .    s 


(3.15) 


or 


TT   (m.  )   _      (m.-l)  (m.  ) 

¥si=si  +A.si 


(W  -   A.I)    s(l)   =  0 
i 

(W  -    A.I)2   s(2)    =   0 


(3.16) 


/  vin.       (m.  ) 

(W-A.IJis      l      =0 


Thus,  s    ,  1  <  j  <m.,  is  a  principal  vector  of  degree  j  associated  with 

,     0 .    0  .       ...     . ,  n     ,      (l )        (m.  ) 

A..   Since  S  is  nonsmgular,  the  m.  principal  vectors  s    ,  ...,  s   l 

are  independent.   If  we  consider  all  the  Jordan  blocks  J   ,  1  <  i  <  £,  we 

m.     - 

l 

have  the  2n  principal  vectors  s  ,  s  , 


,  s^  which  constitute  the 


2n 


similarity  transformation  S  in  (3.12).   Further,  we  rearrange  the  columns 
of  S  such  that  the  first  n  vectors  are  the  principal  vectors  corresponding 
to  the  n  eigenvalues  with  positive  real  parts  of  W. 


10 


Let  us  now  split  J  of  (3.12)  into 


J  = 


:  J2^ 


•  \ 

(3.17) 


where  J     and  J      are   diagonal  with  the  eigenvalues   of  W  with  positive   and 
negative   real  parts   respectively,    and  E     and  E     are  matrices  with  one's  on 
the  super-diagonals    and  zeros   elsewhere.      To   associate  with    (3.6),   let   us 
partition  S  and  S        as 


sll 

S12 

s  = 

S21 

S22 

- 

_ 

s-1- 

Tll 

T12 

T21 

T22 

(3.18) 


where  all  submatrices  are  n  x  n.   Note  that,  by  lemma  2,  it  is  necessary  to 
assume  that  S   and  S       -   S   S     S   are  nonsingular. 


Then 


wt 

e   =  e 

SJS  1t 

=  Se^s'1 

0m 

Sll  S12 

e  x       0 

e^   0 

— 

T 
11 

T 
12 

S21   S22 

n     J2t 
0   e  ^ 

—      — 

0    ^ 

T 
21 

T 
22 

S11eJlteEltT11  +  S12eJ2teE2tT21 


S11eJlteEltT12  +  S12eJ2teE2tT22 


q  eJlteEltT   +  S  eJ2teE2tT     S  eJlteEltT   +  S  eJ2teE2tT 
b21      L   111    22e  d  6  ^  121   b21  ±   X  ^2   b22  d  &   *   X22 


(3.19) 


11 


Since  the  elements  of  Jp  have  negative  real  parts,  e  2  -^  0  as  t  +  »  and 
thus  the  terms  having  e  2  in  (3.19)  vanish  as  t  -*■   °°  .   Therefore, 


X(t) 
Z(t) 


¥t 


T 
C  QC 


„        J_  t  iii_  tm  ^        J.  t  E.,  t_ 

Slie   X  6   1  Tll  Sll6   X  e   2   T12 

21     X  e  X   111  b21  ^   X12 


ri 


T 
C   QC 


Su   eV   eV    (T1X  +  T12   CTQC ) 
S21   eV   eV    (T^  +  T        CTQC) 


(3.20) 


Therefore  the  solution  of  the  matrix  Riccati  equation  is  given  by 

K  =  lim  Z(t)  X_1(t) 
t-*» 

=  lim  [S   eJlt   e¥  (T  +T   CTQC )  ]  [S  eJlW  (T  +T  CTQC)  I"1 

t-K» 


■  S21  Sll 


-1 


(3-21) 


provided  that  the   indicated  inverse   exists.      However,   we   shall   show  that 
only  the  nonsingularity  of  S        and  S        -  S        S  S  p   is  needed  for  the 

stationary  solution  of  the  matrix  Riccati   equation.      From   (3.12), 


T  T 

11        12 


T  T 

21        22 


-1   T 
-A        BR     B 


T  T 

C   QC  A 


Sll      S12 


21        22 


=   J 


(3.22) 


Equating  the    (2,1)   elements   on  each  side  of   (3.22), 


(-T21  A  -  T2£1  CTQC)   Su  +    (T       BR^B1  +  ^  AT)   8^  =  0 


12 
or 

■         -T21   ASn   +   T22   AT  S21   +  T21   BR^  S21+   T22   CTQC   S^   =   0  (3.23) 

Substituting  T21  =  -  T22  S21  S^"1  into    (3.23), 

T22S21Sll'1+   T22ATS21   "   T22S21Sll"lBR"lBTs21   +  ^2°^°  Sll=° 

(3.21+) 


Premultiplying  by  T22"  and  post multiplying  by  S   _1 , 


S21Sll"lA  +  ^^l3!!"1  "  ^l5!!"1^"1^^!3!!"1  +  <?<*   =  °   .0.25) 


which  is  to  be  proved. 
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k.      COMPUTATIONAL  METHODS  AND  NUMERICAL  RESULTS 

The  major  three  steps  for  computing  the  solution  of  the  matrix  Riccati 
equation  are  as  follows: 

1)  Compute  eigenvalues  and  eigenvectors  of  W, 

2)  Compute  principal  vectors  of  ¥  if  W  is  defective, 

3)  Compute  S   S  "   from  S. 

For  computing  the  eigensystem,  we  used  two  methods:   (i)  the  QR 
algorithm  [9]  with  an  inverse  iteration  routine  [10],  and  (ii)  Jacobi-like 
method  [ll]  for  finding  the  eigenvalues  and  eigenvectors.   Their  perfor- 
mances are  satisfactory  and  two  numerical  test  examples  are  given  in  Test 
Result  1  and  Test  Result  3.   When  the  matrix  W  is  nondiagonable,  Test  Result 
2,  the  first  method  requires  the  evaluation  of  the  principal  vectors  as 
well  [12].   This  poses  some  numerical  difficulties  as  we  are  essentially 
trying  to  solve  the  ill-conditioned  system  of  equation 

(W  -  Al)k  x    =  x.,  k  >  1 
where  ^  is  an  approximation  of  true  eigenvalues  A.   Thus,  the  Jacobi-like 
method  is  more  attractive  in  this  case  for  providing  an  approximate  eigen- 
system that  is  suitable  for  engineering  purposes. 

Appendix  I  contains  the  listings  of  the  ALGOL  programs  used  on  the  B-6500 
Appendix  II  contains  the  listings  of  two  ILLIAC  IV  programs,  namely,  the 
inverse  iteration  routine  and  a  routine  for  solving  a  system  of  complex 
linear  equations,  both  written  in  GLYPNIR.   ILLIAC  IV  programs  for  the  QR 
algorithm  [13,  1^]  and  the  Jacobi-like  algorithms  [15,  l6]  are  already 
available. 


Test   Result   1. 


ik 


A  = 


1 

0 

0 

0 

0 

1 

0 

-1 

0 

0 

0 

0 

-1 

0 

0 

0 

0 

1 

0 

-1 

0 

0 

0 

0 

-1 

p  = 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

B  = 


1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

Q  = 


0 

0 

0 

0 

0 

0 

10 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

10 

0 

0 

0 

0 

0 

0 

c  = 


1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

R  = 


1 

0 

0~ 

0 

1 

0 

0 

0 

1 

The  numerical   example  was  tested  by  the  QR  algorithm  and  the   inverse 
iteration  method.      The  ten  distinct   eigenvalues   of  W  are 

±    1. 

±   1.728760U77185  ±  1.577533767573i, 

±  1.35319578U0U6  ±  1.1537U989928Ui. 
We  choose  the  five  eigenvectors  corresponding  to  the  five  eigenvalues 
with  positive  real  parts  and  computing  K  =  S   S     from  them  gives  the 
solution  K  as 
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1.262782609 

2.1+9^009759 

-0.819173651 

0.668267901 

-O.UU3608958 


2.1+9^009759 
7.^  3  5^51161+ 

-1.8257^1858 
1.1229101+32 

-0.668267901 


-0.819173651 
-1.82571+l858 
1.63831+7303 
1. 8257^1858 
-0.819173651 


0.668267901 
1.1229101+32 
1.8257^1858 
7.1+ 35^51161+ 
-2.1+91+009759 


-O.Ul+3608958 
-0.668267901 
-0.819173651 
-2.1+91+009759 
1.262782609 


-1   T 
All  the  eigenvalues   of  G  =  A  -  BR     B  K  have  negative   real  parts   and  hence 

the   above   solution  K  is   acceptable. 


Test   Result   2. 


A  = 


-3       '2 

-2       1 


B  = 


C  = 


1        0 
0        1 


P=Q= 


0        0 
0        0 


The  mat 

rix 

— 

3 

-2 

0 

0 

2 

-1 

0 

1 

¥  = 

0 

0 

-3 

-2 

0 

0 

2 

1 

R  = 


CO 


is   defective,   that   is,   W  has    four  eigenvalues   ±1    (double)   and  each  +1   and 
-1   correspond  to  a  quadratic   elementary  divisor,    respectively.      This 
example  was  tested  both  by  Jacobi-like  method  and  by  a  combination  of  the 
QR     algorithm  and  the   inverse   iteration  method. 
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The  Jacobi-like  method  gives   K  as 
-10.0000081*  6.00000U3 

6.000003^  -k. 0000027 

and  the   QR  and  inverse   iteration  method  give  K  as 
-9. 9999917  5.99999^7 

5.9999953         -3.9999983 
The   exact   solution  is 

~Vlo      6 
6     -k 

For  computing  the  principal  vector  by  the  QR  and  inverse  iteration,  we 
chose  the  perturbation  for  the  true  eigenvalues  as  much  as  10 


Test  Result  3. 


A  = 


6 

k 

k 

1 

h 

6 

1 

k 

h 

1 

6 

k 

1 

h 

h 

6 

B  = 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

C=R= 


1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

P=Q= 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

The  exact  eigenvalues  of  W  are 

±15,  ±5,  ±5,  ±1 
Since  the  matrix  W  is  derogatory,  +5  and  -5  correspond  to  a  two-dimensional 
subspace,  respectively.   The  computed  eigenvalues  by  Jacobi-like  method  are 


IT 


lU. 99999999930 
k.  9999999991+9 

-U. 99999999982 
0.99999999997 


-1^.9999999984 
4. 99999999981 
_1+.  999999999^8. 

-0.99999999995 


The  computed  solution  K  is  given  by 

-1.999999999811  1.9999999999^2  1.999999999913 
1.999999999884  -2.000000000015  -2.000000000000 
2.0000000002T6   -2.000000000^07  -2.000000000378 

-2.00000000036U   2.000000000509  2.000000000I180 


-2.000000000015 
2.000000000087 
2.000000000^80 

-2.000000000568 


in  -which  the  absolute  error  is  less  than  10 


-9 
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APPENDIX  I 


ALGOL  PROGRAMS 
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I    TWO  ALGflL  PROCEDURES'  RlCCATU  AND  RICCATI2.  FOR  COMMUTING  THE  OOOOOlOO? 

I!  SOLUTION  OF  THE  MATRIX  RICCATI  EOUATION  ARE  LISTEO  BtLOWt  FOR  00000200? 

%    "RICCATlU.  THE  PROCEDURE  "EIGENw  (BY  EflERLEIN  AND  bUOTHKoyD)  WAS  00000300? 

%    USEU  WHILE  A  COMBINATION  OF  "HQR"  (BY  MaRTINi  ET  AL)  AND  00000400? 

I  "INVITErATION"  WAS  APPLIED  TO   «RICCATI2%  LINEAR  EqUATIUn  SOLVER  OOOOO5OO? 
%    (BY  BOWDLER*  ET  AL )  WAS  ALSO  USED  TO  PROVIDE  THE  INytRSE  OF  A  MATRIX   00000600? 

%    ANO  TO  CALCULATE  THE  PRINCIPAL  VECTORS  fROM  ( ( W-LAMflUA*! )**K>*X«0.  00000700? 

00000800? 

00000900? 

PROCEDURE  RlCCATIl(N.L»M»A»B#C»P»Q»R«TMX»W»K)l  00001000? 

VALUE  N»L#M,TMXI  00001100? 

INTEGLR  N,L.M.TMX|  00001200? 

ARRAY  A.e(C.P»Q»R»W.K£i»n|  00001300? 

COMMENT  SOLVES  THE  MATRIX  RlccATI  EQUATION  .  OOOOUOO? 

KA*(AT)K-KR(RINV)(BT)K*(CT)QC«o,  00001500? 

A  Is  N*N»  R  IS  N*L»  C  IS  M*N.  P  AND  0  ARE  M*M.  R  IS  L*L  OOOOUOO? 

MATRICES*  RESPECTIVELY*  FORM  2N*2n  MATRIx  W  ANU  COMpUTE  00001700? 

EIGENVALUES  ANU  EIGENVECTORS  OF  W*  WHERE  FOUR  N*N  BLOCKS  00001800? 

Op  w  ARE  C 1 • 1 >«-A»  Cl»2)«8(RlNV>(BT)»  ( 2» 1  )■( C  '  >QC»  AND  00001900? 

(2»2)"(AT).  REARRANGE  EI GEN VEcTOH ( Pr INC  I PAL  VEcT0R>  MATRIX  00002000? 

SO  THAT  THE  pIHST  N  COLUMNS  CORRESPOND  TO  TO  THE  EIGENVALUES  00002100? 

WlTH  P0SITIVE  REAL  rArTS.  dENoTINq  ThE  UPP^R  HaLf  uF  N  COLUMN  00002200? 

VECTORS  BY  SMTXl  AND  THE  LOWER  HALF  BY  SMTX2*  I  HE  SOLUTION  00002300? 

K  IS  THEN  GlvEN  BY  «■( SMTx2) ( SMTxl INv ) .  THE  PRUCEOUrE  00002400? 

"FI6EN"  IS  USED  For  EIGENPrOBlEM  qF  Wi  00002500? 

00002600? 

REGIN  COMMENT  FIRST  OECLARE  PROCEOuRESl  00002700? 

00002800? 

PROCEUURE  SETUP(NtL»M.A»B»C»Q»RlNV.W)|  00002900? 

VALUE  N.L.M,   INTEGER  N.L»M)  00003000? 

REAL  ARRAY  A»B»C»Q.RINV»W[l#l]j  00003100? 

COMMENT   THE  PROCEDURE  SETS  UH  ThE  2N*2N  mATRIX  W,  rInV  Is  INVERSE  OF  RI00003200? 

«EGIN  00003300? 

REAL  ARRAY  RINVBTC1 lL»l IN].QCC1 IM»1 IN] »  00003400? 

REAL  SUMj  00003500? 

INTtGER  I.J.KJ  00003600? 

FOR  Ij"l  STEP  1  UNTIL  N  DO  00003700? 

FOR  J l ■!  STEP  1  UNTIL  N  DO  00003800? 

BtGIN  00003QOO? 

W[I,J]I«-A[  I. J]}  00004000? 

W[N*I,N*J]IbA( J.I1I  00004100? 

ENO)  00004200? 

FOR  1 1 » 1  STEP  1  UNTIL  L  00  00004300? 

FOR  J»«l  STEP  1  UNTIL  N  DO  00004400? 

BEGIN  StJMl.O.Ol  00004500? 

FOR  K,«l  STEP  1  UNTIL  L  DO  00004600? 

S|jM|«SUM*RINVCl»K]*BC  J#KJ1  00004700? 

RINvBTC I.J]l«SUM»  00004800? 

ENOi  00004900? 

FOR  I|M  STEP  1  UNTIL  N  00  00005000? 

FOR  Jl»l  sTEp  1  UNTIL  N  00  00005100? 

BtttIN  SuHlaO.Ol  00005200? 

FOR  Ki«l  STEP  1  UNTIL  L  DO  00005300? 

SuH«"SUH*B[IsK]*RINVBTtK#J]|  00005400? 

Wt I,N*J]l«SUMl  00005500? 

END|  00005600? 

FOR  I|»l  STEP  I  UNTIL  M  DO  00005700? 
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FOR  Jl«l  STEP  1  UNTIL  N  DO 

BEGIN  SuMl.O.Oi 

for  K|«i  step  i  until  m  do 

SuMl"SUH*OtItK]*CCK.JJI 
QCC I #  J) I -SUM l 
ENO| 
TOR  I|«l  STEP  1  UNTIL  N  DO 
FOR  J|«l  STEP  1  UNTIL  N  DO 
BtG]N  SuNlsO.Ol 

FOR  K|»l  STEP  1  UNTIL  N  DO 
StjMI»SUM*CtK.  i  J*«C  t  K.  J3 1 
W[N*I.J]|"SUMI 
ENDI 
ENO  OF  SETUP| 

COMMENT  INSERT  THE  PROCEDURES  INNPRODtDECOMPOSE.  AND  ACCSOLyEi 

PROCEUURE  lNVERSE(N»AtX)l 

VALUE  N|  INTEGER  Nj 

DOUBLE  ARRAY  A»XC  1  • 1 1  I 

COMMENT  THE  PROCEDURE  COMPUTES  THE  INVERSE  OF  AN  N*N  rUl 

UNSyMMcTRIc  MATRIX  A.  USINq  ThE  PROCEOURgS  "DECQMpUsE*  AND 
"ABSOLVE".  WRITING  ON  X| 
BEGIN 

INTtGER  I.J.IT.D2I 
REAL  Oil 

ARRAY  AA'BB.BU I N* 1 iNJ.P[llNll 
FUR  I • * 1  STEP  1  UNTIL  N  DO 
1  UNTIL  N  Oo 

BCl»J]l«ltO  ELSE  B(I»J]l«O.OI 
1  UNTIL  N  DO 
1  UNTIL  N  DO 


FuR  Jl«l 

IF  J«I 

FUR  I,«l 

FUR  J  t -1 


STEP 
THEN 
STEP 
STEP 
AACT.Jlt-ACI.J1l 
DtCOMpOSE(N»AA.Dl.02»P)l 
AtCSULVE(N.N.A,AA.P»B.X»BB»IT)| 
ENO  INVERsEl 

COMMENT  iNjSfRT  THE  PRnCEDURE  EIGENi 

PROCEUURE  MuLT(Nl.N2iNJ»A.B»C)l 
VALUE  N1»n2.N3I  INTEgEH  N1*N2»n3> 
DOUBLt  ARRAY  A.B.Ctl»lJI 
REGIN 

INTtGER  I.J.KI 

REAL  SUMI 

FOR  I l » 1  STEP  i 

FOR  J  t • 1  STEP  1 

BEGIN  SUM|«0.0| 

FUR  K,"t  STEP 

SUMI»SUM*ACI»K]*BCK,J]| 
CI  I  .J]  l«SUMI 
ENO» 
FNO  MULTl 

ARRAY  HiNv»SMTXl,SMTX2lllN.llN].WW.TEMP,VECT[llN*N»llN*NJ| 
INTLGEH  ArRAy  search(1»nii 
INTtGER  N?» I. J- INDEX' 


UNTIL 
UNTIL 


Nt 
N3 


DO 
OO 


1  UNTIL  N2  OO 


00005800? 
00005900? 
00006000? 
00006100? 
0000*200? 
00006300? 
OOOO64OO? 
00006500? 
00006600? 
00006700? 
00006800? 
00006900? 
00007000? 
00007100? 
00007200? 
00007300? 
00007400? 
00007500? 
00007600? 
00007700? 
00007800? 

00007900? 
00008000? 
00008100? 
00008200? 
00008300? 
00008400? 
00008500? 
00008600? 
00008700? 
00008800? 
00008900? 
00009000? 
00009100? 
00009200? 
00009300? 
00009400? 
00009500? 

00009600? 
00009700? 
00009800? 
00009900? 

00010000? 
00010100? 
00010200? 

00010300? 

00010400? 
00010500? 

00010600? 
00010700? 
00010800? 
00010900? 

00011000? 
oooiuoo? 

00011200? 
00011300? 
00011400? 
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N2I«N+N|  00011500? 

INVtRSE(L.R»RlNV)l  0001U00? 

SETUP <N,L#M»A»B»C»Q»RINV»W)J  00011700? 

TOR  It*l  STEP  1  UNTIL  N2  DO  00011800? 

FOR  Jl "1  STEP  1  UNTIL  N2  DO  00011900? 

TEMPI  I. JJ  t«M[ I.J] I  00012000? 

EIQtN(N2»TMX,TEMP,VECT)l  00012100? 

COMMENT  SEARCH  N  COLUMNS  OF  VECT  CORRESPONDING  TO  N  EIGENVALUES  00012200? 

WITH  POSITIVE  REAL  PARTS.  SMTXl  AND  SMTX2  ARi  IT*  UPPER  00012300? 

AND  LOWER  HALVES*  RESPECTI VELY|  00012400? 

INDLX|.0l  00012500? 

FOR  I » » 1  sTEp  1  UNTIL  N2  DO  00012*00? 

I*     TEMP[I»I1  GTR  0  THEN  00012700? 

BLGIN  00012800? 

INDEX|«INDEX*1 ,  00012900? 

SEArChI INDEX] tail  00013000? 

ENO|  00013100? 

FOR  J|»1  STEP  1  UNTIL  N  DO  00013200? 

FOR  I, al  STEP  1  UNTIL  N  DO  00013300? 

BEGIN  00013400? 

SMTxltliJ) |«VECTC I » SEARCH! J]]j  00013500? 

SMTX2tI.J]l«VEcTtN*I»SEARcHCJ}]l  00013600? 

ENDI  00013^00? 

INVLRSE(N, SMTXl. TEMP)I  00013800? 

MULHN»N»N»SMTX2.TEMP»K)|  00013900? 

END  RiCCATH»  00014000? 

00014100? 

00014200? 

PR0CEUURE  RICCatI2(n.L»M»a»B»C»P»0»R»W»K»0EFECt»DER0(»AIE)I  00014  300? 

VALUE  N.L.Mf  00014400? 

INTEGtR  N.L.Ml  00014500? 

ARRAY  A»3.C.PfQ»R*H*Ktl'l)l  00014600? 

INTEGtH  ARRAY  DEFECT, DEKUGATEt 1 J  J  00014700? 

COMMENT  SnLvES  THE  MATRIX  RICCATI  EQUATION  00014800? 

KA*(AT)K-KB(HINV)(BT)K*(CT)QC«0.  00014900? 

A  Is  N*N.  b  IS  N*L»  C  IS  M*N.  P  AND  Q  ARE  M*M,  R  IS  L*L  00015000? 

matrices*  respectively.  form  2n*2n  matrix  w  and  compute  00015100? 

eigenvalues  anu  eigenvectors  of  w,  where  four  n*n  blocks  00015200? 

OF  W  ARE  (l.l)«-A.  (1.2)«B(RINV)(BT)*  (2*1)«(CU0C  AND  00015300? 

(2»2)«(AT).  REARRANGE  El GEnVEcTOR( PRlNC I  PAL  VECTOR*  MATRIX  00015400? 

SO  THAT  THE  FIRST  N  COLUMNS  CORRESPOND  TO  TO  THE  EIGENVALUES  00015500? 
WITH  POSITIVE  HEAL  PARTS.  DENOTING  THE  UPPER  H*LF  Up  N  COLUMN    00015600? 

VECTORS  BY  SMTXl  AND  THE  LOWER  HALF  BY  SMTX2,  THE  SOLUTION  00015700? 

K  IS  THEN  GIVEN  BY  K«( SMTX2) < SMTXl INV ) .  THE  PROCEDURES  00015800? 

"hShLD"*  "hOR**  AnD  "INVITEKATION"  ARE  USEq  FOK  EH»ENPRORlEM  00015900? 

OF  Wl  00016000? 

00016100? 

REGIN  COMMENT  FIRST  DECLARE  PROCEDURES!  00016200? 

00016300? 

PROCEUURE  SETUP(N»L*MtA*B,C»Q*RINV*W)|  OOOI64OO? 

VALUE  N.L.Mf   INTEGER  N.LtM*  00016500? 

REAL  ARRAY  A  »  B .  C  » Q .  RIN V  •  W  [  1  i'l  )  I  00016600? 

COMMENT  THE  PROCEDURE  SETS  UP  ThE  2N*2N  MATRIX  W.  RINV  IS  INVERSE  Of  R!  00016700? 

QEGIN  00016800? 

REAL  ARRAY  R I NVBT[ 1  I L« 1  I N J . QC [ 1  I M» 1 1 N ] I  00016900? 

REAL  SUM1  00017000? 

INTLGER  I.J.KI  00017100? 
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FOR  I t -1  STEP  1  UNTIL  N  DO 
TOR  J|»1  STEP  1  UNTIL  N  DO 
BtGIN 

M[I.J]I«-ACI»J]| 

WtN+I.N*J]l«ACJ#ni 

ENDI 

Fo«  i»«i  step  i  until  l  Do 

FOR  Jial  STEP  1  UNTIL  N  DO 
BtGIN  SuM|aO»0| 

FOR  K|«l  STEp  1  UNTIL  L  DO 
SuM|«SUM*RJNVt I»KJ*B[ J»K]  J 

RINvBTtI»J]l«SUMI 

END| 

FOR  I|»l  STEP  1  UNTIL  N  00 
FOR  J  t"l  STEp  1  UNTIL  N  DO 
BtGIN  SuMl-0,01 

FOH  Kt«l  STEP  1  UNTIL  L  00 

SuM,«SUM*BU.K)*RlNVBT[K.JJ| 
W[I.N*JJ|«SUmI 
ENO| 
FOR  Iial  STEP  l  UNTIL  M  00 
FOR  J|"l  STEP  1  UNTIL  N  DO 
BLqIN  SUMI.O,0| 

FOR  Kf.l  STEP  1  UNTIL  M  DO 

SUM|.SUM*Q[ I.K]*C[K»Jj I 
«CC  If j] I»SUMJ 
ENDl 
FOR  Il»l  STEP  1  UNTIL  N  DO 
FOR  J  t "1  STEP  1  UNTIL  N  DO 
BtGIN  SUMl«O.0l 

FOR  K|b1  STEP  I  UNTIL  M  DO 
SuM|«SUM*C(K.n*UC[K»Jj; 
W[N*I,j) l«SUMI 
END> 

end  of  setupi 

comment  Insert  the  procedures  innprod. decompose*  and  accsolvei 

proceuure  inverse(n»a#x)| 

VALUE  N|  INtEgEH  N| 

nOUBLt  ARRAY  A.XU.lJI 

COMMENT  ThE  PROCEOURE  COMPUTES  THE  INVERSE  OF  AN  N*N  RtAL 

unsymmetric  hatrix  a.  using  The  procedures  "DECQmpu$e" 
waCcsolvE".  writing  on  XI 

BEGIN 

INTtGER  I.J.IT.D2I 
REAL  nl  | 

ARRAY  AA'BB.BUtN.l  fNJ»Ptl  tNJ| 
FUR  I i - 1  STEP  1  UNTIL  N  DO 

DO 

B[I»J]nO«OI 


AND 


FuR  J  I -l 

IF  J-I 

FUR  I,-! 

FUR  J i ■ 1 


END 


step 

then 
step 

STEP 

AAt !ij]|iA[I,jii 
ntCnMp0sE(N.AA,0l.02#P)| 
ACCSULVE(N#N.A.AA.P»B#X»BB»IT)| 


i  until  n 

fl[I. J]  i-l .0  ELSE 
1  UNTIL  N  DO 
1  UNTIL  N 


00 


00017200? 

00017300? 

00017400? 

00017500? 

00017600? 

00017700? 

00017*00? 

00017900? 

00018000? 

00018100? 

00018200? 

00018300? 

00018400? 

00018500? 

00018600? 

00018700? 

00018800? 

00018900? 

00019000? 

000i9i00? 

00019200? 

00019300? 

00019400? 

00019«00? 

00019600? 

00019700? 

00019800? 

00019900? 

00020000? 

00020100? 

00020200? 

00020300? 

00020400? 

00020500? 

00020600? 

00020700? 

00020800? 

00020900? 

00021000? 

00021100? 

00021200? 

00021300? 

00021400? 

00021500? 

00021600? 

00021700? 

00021800? 

00021900? 

00022000? 

00022100? 

00022200? 

00022300? 

00022400? 

00022500? 

00022600? 

00022700? 

00022800? 
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00022900? 

Comment  InseRt  the  procedures  hshld  and  hori  ooo23ooo? 

00023100? 

PROCEDURE  MuLT(Nl»N2*Ni'A»B«C)l  00023200? 

VALUE  N1»N2*N3|  INTEGER  N1.N2.N3j  00023300? 

DOUBLE  ARRAy  A.B.CU.Ul  00023400? 

REGIN  00023500? 

!NTtOER  I»J»K»  00023600? 

REAL  SUMS  00023700? 

FOR  Il«l  STEP  1  UNTIL  Nl  DO  00023800? 

FOR  Jl«l  STEp  1  UNTIL  N3  00  00023900? 

BEGIN  SUM|«0,0I  00024000? 

FUR  K|«l  STEP  1  UNTIL  N2  DO  00024100? 

SUM i * SUM* A [ I»K]*B[K«J)|  00024200? 

CI  I • Jl l.SUMI  00024300? 

ENOI  00024400? 

END  MULTl  00024500? 

00024600? 

PROCEUURE  lNVlTERATION<N»w»X»KK»P>|  00024700? 

VALUE  N,  TNTEGER  N.KKj  00024800? 

DOUBLE  ARRAy  W.Xtl.lli  00024900? 

ARRAy  P[l]l  00025000? 

CnMMENT  SnLVES   WX.O  BY  INVERSE  ITERATION*  WHERE  W  IS  AN  N*N  REAL  00025100? 

MXTRlX.  WRITING  SOLUTION  ON  X.  THUS.  THE  PROCEDURE  COMPUTES  00025200? 

ThE  EIqENvECTOKS  CORRESPONDING  TO  THE  GIVEN  EIGENVALUES.         00025300? 

N  Is  THE  ORDER  OF  M  AND  KK  IS  THE  NUMBER  OF  iTtRATiONS           00025400? 

NEEqEO  FOR  ITERATING  EACH  EIGENVECTOR!  00025500? 

BEGIN  00025600? 

INTfcGER  I.J»D2»R.L»MAXIN0EX|  00025700? 

DOUBLE  DUNAXVALi  00025800? 

DOUBLE  ARRAY  MN( 1 1 N. i | NJ »B.B8[ i | N» 1 | U |  00025900? 

LABEL  INVIT.OUTI  00026000? 

FOR  Il"l  STEP  1  UNTIL  N  DO  00026100? 

FOR  Jl ■ l  STEP  1  UNTIL  N  DO  00026200? 

ww [  I »  J  ]  t»ri t  I » J  3 1  00026300? 

DECUMpOsE(N»WW»Dl»02»P)l  00026400? 

FOR  Il«l  STEP  1  UNTIL  N  DO  00026500? 

BII,1]»81,0»  00026600? 

KK I "1 1   R| ■ I j  00026700? 

INVITI  00026800? 

ACCSOLVe(N»R»W»WW»P»B»X.BB*L)I  00026900? 

MAXVALUOj  00027000? 

FOR  Ii"l  STEP  1  UNTIL  N  DO  00027100? 

If  ABs(X(I*l]>  GTR  MAXVAL  THEN  00027200? 

BEGIN  00027300? 

MAXvALl«ABS(Xr I#i])l  00027400? 

MaXiNDEXI«P  00027500? 

ENDi  00027600? 

MAXvA|_l«X[MAXlNDEX»ni  00027700? 

FOR  I|«l  STEP  1  UNTIL  N  DO  00027800? 

XlI»l]l«X[Iil]/MAXVALI  00027900? 

FOR  III}  STEP  1  UNTIL  N  00  00028Q00? 

I*     ABS(BCI»1]-Xtl»l])  GTR  P"12  THEN  00028100? 

BLGIN  00028200? 

IF  KK  GTR  10  THEN  GO  TO  OUTj  00028300? 

FOR  J|-l  STEP  1  UNTIL  N  DU  00028400? 

BtJ,l]l"X[j.i ]l  00028500? 
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KK»»KK+U 
GO  TO  INVITI 
END) 

nuTi 

ENO  INVREsElTi 

INTtGER  I.J»JJ»CNT.N2»INDEX#ITC0UNT| 

ARRAY  TEMpfORGNtSHTX*veCTORtt|N«N«l|N4N]*RlNViSHTXUSHTX2tl|N«l|N]» 

EIGR»EIGI»!NTCH(1IN4N)I 
REAL  PERTURB. MPUiERl 
INTtGER  ARRAY  SEARCHllIN)! 
N2l"N*N| 

INVk.RSE(L.R.RlNV)l 
SETUP(N,L,M,A»B,C.Q.HINV.H)| 
TOR  I|«l  STEP  1  UNTIL  N2  00 
FOR  J|"l  STEP  1  UNTIL  N2  00 


HSHLO(i,N,,TEMP)| 
HQR(N2»TEMP.I 


N2 
N2 


DO 
00 


N2  00 


H 


EIQR.EIGI* CNT)| 
FOR  I | »1  STEP  1  UNTIL  N2  00 
FOR  J I  - 1  sTEp  1  UNTIL  N2  00 

ONGHCli J]l«NCItJ]l 
COMMENT  COMPUTE  EIGENVECTORS  OR  PRINCIPAL  VECTORSi 
FOR  Jjl-1  STEP  1  UNTIL  N2  00 
BEGIN 

FUR  I t»l  STEP  l  UNTIL 

FUR  Ji-i  STEP  1  UNTIL 
*[I»J] ««ORGW[I,J]l 

FUR  I » - 1  STEP  1  UNTIL 

•»tIiI]l"WU.i]-ElGR[  JJ]*PERTURB| 

JJ«DEFECT[JJ]  then 
BEGIN 

FOR  Ii-1  STEP  1  UNTIL  N2  00 

W[I.I  J»«Wti,IJ*pERTURB*MpLIERl 

MuLTCN2.N2»N2»W»W,TEMP)l 

Ff)R  Ital  STEP  1  UNTIL  N2  00 
FOR  Jlal  STEP  1  UNTIL  N2  DO 
*t  I«J]i«TEMP[  I#J)I 
END  | 
IF  JJaOfROGATEUJ]  THEN 

FOH  Ii-l  STEP  1  UNTIL  N2  00 
Mf I • 1 1 |»WC 1. 1 1 ♦PERTURB* J Jl 

INVITFRATI0N(N2#W.VECT0R»ITC0UNT»INTCH)» 
FUH  I  t  a.j  STEP  i     UNTIL  N2  DO 

FUR  Jiaj.2  DO 

SMTx(l»JJ]|avEcTORtI#l Ji 
END  JJ> 

COMMENT  SEARCH  N  COLUMNS  OF  SMTx  CORReSpONOiNG  TO 
WITH  POSITIVE  REAL  PARTS,  SMTxl  AND  SMTx2 
AND  LOWER  HALVES*  RESPECTI VELYl 
INDtXiaQi 

FOR  ll«l  STEP  1  UNTIL  N2  00 
l>     ElGRfl]  GTR  0  THEN 
BLGIN 

iNOFXtalNUEX^l  | 
SEArcmC INDEX]  I.I  > 
ENDl 


N  tlGENVALUES 
ARt  ITS  UPPER 


00028600? 
00028700? 
00028800? 
00028900? 
00029000? 
00029100? 
00029200? 
00029300? 
00029400? 
00029500? 
000296O0? 
00029700? 
00029800? 

00029900? 
00030000? 
00030100? 
00030200? 
OOO3O1OO? 
00030400? 
00030500? 
00030600? 
00030/00? 
00030800? 
00030900? 
00031000? 
OOO3UOO? 
00031200? 
00031300? 
00031400? 
00031500? 
00031600? 
00031700? 
00031800? 
00031900? 
00032000? 
00032100? 
00032200? 
00032300? 
00032400? 
00032500? 
00032.00? 
00032700? 
00032800? 
00032900? 
00033000? 
00033100? 
00033200? 
00033300? 
00033400? 
00033500? 
00033600? 
00033700? 
00033^00? 
00033900? 
00034000? 
00034100? 
00034200? 
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FOR  Jl«i  STEP  1  UNTIL  N  DQ  00034300? 

FOR  Ii"l  STEP  1  UNTIL  N  DO  00034400? 

BEGIN  00034500? 

SMTXltI.J]l«SMTX[I»SEARCH[Jlll  00034600? 

SMTX2(I,J)I«SMTX(N+I'SEARCH(J])J  00034700? 

CNDI  00034800? 

INVtRSE(N,SMTXl»TEMP)l  00034900? 

MULT(N»N»N»SMTX2.TEMP»K)|  00035000? 

END  RlCCATl2'  00035100? 
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APPENDIX  II 


ILLIAC  IV  PROGRAMS 


29 


*  TWO  ILLIAC  IV  GLYPNIR  PROGRAMS.  "ClINSYs"  AND  "INVITE-RATION"  ARE 

f  LISTED  BELOW,  "CLlNSTS"  IS  USED  TO  OBTAIN  THE  INVERS*  OF  A  MATRIX 

%    AND  ThE  SOLUTION  Or  A  COMPLEX  SYSTEM  OF  UNSYMMETRIC  LINEAR  EQUATIONS. 

%    MNVITERATION"  IS  USED  TO  CALCULATE  THE  PRINCIPAL  VECTORS  OF  DEGREE 

*  K  FKOM  ((WLAMBDAM)**K)*XnO. 

SUBROUTINE  cLINSyS(CINT  N.ClNT  R#PCPOINT  A'PCPOINT  x»PCPOINt  B)J 
BEGIN  x  THE  SUBROUTINE  SOLVES  A  COMPLEX  SySTEM  OF  UNSyMMETKlC  LINEAR 

x  Equations  ax»r.  a  is  n*2n  and  b  is  n*2R  matrices*  respectively. 

*  ThEir  qEnErAi  EiEMENTS  ArE  At  I  »2J-H*I*  At  I  »2J]  AND  BCp2j-ll* 

*  I*B[I,2J],  RESPECTIVELY,  R  IS  THE  NUMBER  OF  RIGHT-HAND  SIDES, 

SUBROUTINE  INNpROOf  CREAL  CWCREAl  C2»PREAL  AK»PREAL  BK»C«EAL  Dl# 

CREAL  D2)» 

begin  *  accumulatEs  the  sum  of  products  and  aods  it  10  the  initial 
%   value  (c1.c2j. 

CHEAL    Sl2| 

S12,»CUC2| 

Si2l»Sl2*ROwSUM(AK*BK)i 

01  I .Si  2  t 

D2I.S12-D1I 

END) 

SUBKgUTTNE    CDEC0MP0SE(CiNT    N»PCPnlNT    A'CREAL    DErR'CRtAL    i>etI» 


CI 
BEGAN    %    ThE    COMPLEX    UNSy 

where  l  is  lower 

MATRICES.  AND  Gv 
AND  ITS  GENERAL 
KEEpS  THE  RECOSO 
A,  THE  OETERMIMA 

C1NT  I» j.K.L.P.PPj 

PREAL  VICTOR  ATI6AJI 

PKEAL  Z7I 

CKEAL  X.Y.Z.V.W.H.HHI 

LABEL  FAlLl 
MUDEI«TRUEl 

MUDE»«TRUE  AND  pen  lss 
LUUP  II»1.1»N  DO 

lNNpRoD(0»O.AtI]»AtI 
DLTRI.ii  DETII.Oj  DETE 
LUOP  Kl-l.l.N  00 
BLGIN 

Lf«Kl  Pl«K*K|  PPI«P- 

LOOP  II«K»1.N  DO 
BEGIN 

MoDEI«TRUEl 

MoDE»«TRUE  AND  PEN 

Innprodc-graboneca 
InnpRod(-h,-hh.rtl 
Innprodcgraroneu 

lNNpROD(H.HH»ACIl. 
Y|«-HJ 

modei«true» 
m0de1.true  and  pen 
mode»«truej 


NT  DETE»CNP0INT  INT)  I 

MMETRIC  MATRIX  A  IS  DECOMPOStD  INJO  LU. 

triangular  and  u  is  unit  upper  triangular 
erwritten  on  a,  a  is  stored  in  N*2N  array 

ELEMENTS  ARE  At I*2j«l 1«I*A( I *2 J] t  INTIM 
Of  ANy  INTERCHANGES  mAdE  To  THE  ROWS  Op 
NT  (0ETR*I#0ETI)*2**DETE  IS  ALSO  COMPUTFD, 


N  +  NJ 

]ilNT(I].H)l 
1*0 1 


1)  Zl»OJ 


LSS  N*N| 
tI].PP"i)*0*A[I]»ATtPP].H*HH^ 
(l»*A(l]).ATtP].X,HH)| 
tIJ.P-t),O.RTL(l»»ACI]).ATtPPj.H»HH)l 

AT(Pl.H»HH)l 


EQL  PP"ll   AC  1 3  f «X| 


00000100? 
00000200? 
00000300? 
OOOOOaOO? 
00000500? 
00000600? 
00000700? 
OOOOOBOO? 
00000900? 
00001000? 
00001100? 
00001200? 
00001300? 

oooouoo? 

00001500? 

oooouoo? 
00001700? 

00001800? 
00001900? 
00002000? 
00002100? 
00002200? 
00002300? 
00002400? 
00002500? 
00002600? 
00002700? 
00002800? 
00002900? 
0000300^? 
OC003100? 
00003200? 
00003300? 
00003400? 
00003500? 
00003f00? 
00003700? 
00003800? 
00003900? 
0000*000? 
00004100? 
00004200? 
00004300? 
00004400? 
00004500? 

00004600? 
00004700? 
00004S00? 
00004900? 
00005000? 
00005100? 
00005200? 
00005300? 
00005400? 
00005500? 
00005600? 
00005700? 
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00006800? 
00006900? 
0000?OOP? 


So°oDI!:TTRRUULAN°    PEN    E°L    P'U  AI»-TI                                                                       00005800, 

HODE„TRUE    AND    PEN.II  SSSSfSSSJ 

l'    *    GTR    Z    ™EN  00006200? 

ENB0riN    Z,"X'    L""    ENDl  oZZlll 

lFDL,SKUK,T!SSCUTm    AN°    ^  LSS    "*"*                                                                       OOOoJjSSJ 

BEGIN    DETHI-OFTRI  XXSXJJaXI 

22  j  »ACK3 1 
AfKll-Afi 1» 

iNTtLi'^TNTrKii 

E               *WT|Uj lalNiCKjl  00007200? 

lNTfKli«il  00007300? 

x.-JrJboUcackj.pp-u,  SSJSfJSS? 

li:«™*«t*un»  Z°o7rl°o°ol 

OETj,;JjDETI*V*OETRl  SSSS^SSS? 

IF  ABS(DETR)  GTR  ABS(OETI)  SSSSl?So? 

Then  whdetr  else  hi-detii  §§§§si§S? 

JErG?N"L°ARTEHLNL?EGlN  DETE'"01  G°  T°  FAIU  EN°'                        0S0O83OO? 

Lli    IF  ABSCW,  GEO  1  THEN  22SS2S222- 

BEGIN  W!.W*0.0625,  000860 

°lVATrTi:°a\06'5t  SSSS55SSJ 

DETE.-DETEMI  00008800? 

EN0|   G°  T°  L1»  00008900? 

EN0  00009000? 

BEGIN  LABEL  L2i  00009100? 

L2,    IF  IbjCO  LSS  0.0625  THEN  SSSSSSSS5 

BEGIN  S^ETR-H,  =  °0' 

DETKI-0ETR*l6l  00009500? 

0ETII.0ET1M6I  00009*00? 

GO  TO  l\]  00009700? 

end, 

MnOEl.TRUE  ANO  PEN  LSS  K*K-2|  00010*00* 

lNNPROO(-GRABONE(AtK).PP-n.O.AtKl.ATtPP],H»HH)l  00010600? 

INNP«nO(-H.-HH.RTL<J..A[K]).AT[Pl,V.MH)l  00010700? 

|NlJpR0D(-6RAHnNE(ACK].P-n,0.RTL(l,.A|RPl).AT[PPl,H.MH)l  OOOIOSOO? 

INNPR00(H.HHPA[K3.AT[P3.H.HH)|  OOOloSooJ 

SaScVimucf  SSSInSS? 

MROEI-TRUE  ANO  PEN-PP-U  OOOllJooJ 

Mnon'THun  oooiuoo? 


0000980O? 
00009900? 
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modemtrue  ANO  »EN«P-1I  00011500? 

AfK]|»<W*X-V*Y>/Zl  00011600? 

END|  00011700? 

ENO|  00011800? 

FAIL!  00011900? 

ENOl  «0r  CDECOMPOSE.  00012000? 

00012100? 
SUBHOUTINE  CACCS0LVE(CINT  N.CINT  R»PCP0INT  A.PCPOINT  AA.CNPOINT  P»     00012200? 

PCPOINT  B.PCPOINT  X.PCPOINT  BB.CiNT  I) |  00012300? 

BEqIN  *  SOLVES  AX«B.  WhERE  A  IS  AN  N*2N  cOHPlEX  UNSY"NETKIC  ANq  00012400? 

»  b  is  an  n*?r  complex  matrices*  respectively.  aa  1$  lu  00012500? 

*  decomposition  proouceo  by  "cdecompose".  the  hesiuuals  00012600? 
%   bb«b-ax  are  also  calculated  and  ad*88  is  solved.  over-      00012700? 

*  Writing  d  on  bb.  l  is  the  number  of  iterations.  00012800? 

00012900? 

SUBROUTINE  CS0LVE(CINT  N.CINT  R.PCPOINT  A.CNPOINT  iNT.  00013000? 

PCPOINT  B)|  00013100? 

BtGIN  00013200? 

CINT  I,J,K,KK,P,PPl  OOOI33OO? 

PREAL  VECTOR  BT[6«]I  00013400? 

CREaL  X.Y»Z.21.Z2»H.HH»  00013500? 

PREAL  VECTOR  Ct6«l|  00013600? 

BOOLEAN  AMOOEI  00013700? 

LouP  I««1»1.N  On  00013800? 

IF  INTCII  NEQ  I  THEN  00013900? 

LOOP  Jt«R*R."l»l  00  00014000? 

begin  oooiaioo? 

Xi«QRABONE(B£IJ»J-l)l  000U200? 

MODEI-TRUE  AND  PEN"J-ll  00014300? 

B[Ilt»GRABONE<BClNT[U].J-l)l  0001 4400? 

Bt  INTC 111 l"XI  00014500? 

EnOI  00014600? 

LOOP  K»«R*R»-2.2  DO  00014700? 

BEGIN  00014800? 

KkIbK-11  00014900? 

LOOP  Il'l.l.N  00  00015000? 

BEGIN  00015100? 

MODEl'TRUE  ANO  PEN  LSS  1*1-2 J  00015200? 

lNNPROO(-GRABONE(B(I]*KK-l).0»AtI]»BT(KK)»H.HH)|  00015300? 

lNNPR0DCH,-HH.RTL(l..AHl).BTtK].X.HH)l  00015400? 

INNPR0D(-GRAB0NECBII1»K-1)»0#RTL(1»»A[I]).BTIKK]»H»HH)J  00015500? 

INNPR00(H.HH»A[I J.BTtKl.H.HH))  00015600? 

Y|«-HI  00015700? 

P|«I«II  PP|«P-i|  00015800? 

Zi!»GRABONE(  ACIJ.PP-DI  00015900? 

Z2»»GRABoNF(AU J.p-1 >1  00016000? 

ZiaZl*Zl«Z?*42l  00016100? 

MoOEl«TRUE  AND  PEN.KKJ  00016200? 

Btl  Jl.(X#Zl*Y*Z2)/Zl  00016300? 
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